
%% Compute Steady-State

in = 0 ; % Initial infected fraction
d  = 0 ; % Initial dead fraction 
s  = 1 ; % Initial susceptible fraction
r  = 0 ; % Initial recovered fraction

delta = deltascale*(deltabar + exp(phi * (kappa * in))-1) ;

for j=1:num_sectors
eval(['m', int2str(j), ' = 0 ;']) ;
eval(['M', int2str(j), ' = 0 ;']) ;
eval(['l', int2str(j), ' = 1 ;']) ;
eval(['L', int2str(j), ' = (1-d-kappa*in)*l', int2str(j), ' ;']) ;
eval(['mbar', int2str(j), ' = 0 ;']) ;
eval(['chi', int2str(j), 't = chibar', ...
    int2str(j), '*(1-Deltachi', int2str(j), '*(1-exp(-rhom*mbar', int2str(j), '))) ;']) ;
eval(['M', int2str(j), ' = 0 ;']) ;
eval(['c', int2str(j), ' = (1-d-kappa*in)*(l', ...
    int2str(j), ' - chi', int2str(j), 't/2*m', int2str(j), '^2) ;']) ;
eval(['C', int2str(j), ' = (N-d)*c', int2str(j), ' ;']) ;
eval(['P', int2str(j), ' = z', int2str(j), '*(c', int2str(j), '/C', int2str(j), ')^(-1/sigma) ;']) ;
end

c = 0 ;
l = 0 ;
for j=1:num_sectors
    eval(['c = c + (1/num_sectors)*z', int2str(j), '*c', int2str(j), '^((sigma-1)/sigma) ;']) ;
    eval(['l = l + (1/num_sectors)*(1/(1+eta))*l', int2str(j), '^(1+eta) ;']) ;
end
c = c^(sigma/(sigma-1)) ;

e_by_c = 0 ;
e_by_l = 0 ;

for j=1:num_sectors
    eval(['e_by_c = e_by_c + (1/num_sectors)*(1-d)*ec', int2str(j), ...
        '*c', int2str(j), '*C', int2str(j), ';']) ;
    eval(['e_by_l = e_by_l + (1/num_sectors)*(1-d-kappa*in)*(1-m', int2str(j), ...
        ')*el', int2str(j), '*l', int2str(j), ...
        '*L', int2str(j), '*(1-M', int2str(j),') ;']) ;
end

e       = e0 + e_by_c + e_by_l ;
lambdae = ((1/c)-l1^eta)/(el1*(1-m1)*(1-M1)*L1+ec1*C1) ;
lambda  = (1/P1)*(l1^eta + lambdae*el1*(1-m1)*(1-M1)*L1) ;

lambdadtmp = 0 ;
for j=1:num_sectors
eval(['lambdadtmp = lambdadtmp - (1/num_sectors)*lambda*(l', int2str(j), ...
    '-chi', int2str(j), 't/2*m', int2str(j), '^2 - c', int2str(j), ...
    ') + lambdae*(ec', int2str(j), '*c', int2str(j), '*C', int2str(j), ...
    '+ (1-m', int2str(j), ')*el', int2str(j), '*l', int2str(j), '*L', int2str(j), '*(1-M', int2str(j), ')) ;']) ;
end

lambdad = (1/((1-(1/beta)))) * ( l - ud - log(c) + lambdadtmp) ;%(C^(1-sigma)-1)/(1-sigma)

lambdaitmp = 0 ;
for j=1:num_sectors
eval(['lambdaitmp = lambdaitmp -  (1/num_sectors)*kappa*lambda*(l', int2str(j), ...
    '-chi', int2str(j), 't/2*m', int2str(j), ...
    '^2) + kappa*lambdae*((1-m', int2str(j), ')*el', int2str(j), '*l', int2str(j), ...
    '*L', int2str(j), '*(1-M', int2str(j), ')) ;']) ;
end

lambdai = (1/((1-rho-delta*kappa+R0*rho*s-1/beta))) * ( kappa*l - us*kappa  + lambdaitmp - delta*kappa*lambdad ) ;
lambdas = 0 ;
    
Vmtmp = 0 ;
for j=1:num_sectors
eval(['Vmtmp = Vmtmp + P', int2str(j), ...
    '*chibar', int2str(j), '*Deltachi', int2str(j), ...
    '*exp(-mbar', int2str(j), ')*m', int2str(j), '^2 ;']) ;
end
    
for j=1:num_sectors
eval(['Vmbar', int2str(j), ' = (1/(1-beta*psim))*(1/2)*(1/num_sectors)*lambda*(1-d-kappa*in)*Vmtmp ;']) ;
end

Vi = -(1/beta)*lambdai ;
Vs = -(1/beta)*lambdas ;
Vd = -(1/beta)*lambdad ;

gamma    = R0*rho/e ; % Transmission rate for given R0

params  = [ ...
	beta	 
	sigma	 
	N		 
	gamma	 
	deltabar 
	phi		 
	kappa	 
	rho		 
	e0		 
	eta		  
	us		  
	ud		  
	psim 	  ] ;

for j=1:num_sectors
params = [params ;
	eval(['ec', int2str(j)]) ;
	eval(['el', int2str(j)]) ;
	eval(['chibar', int2str(j)]) ;
	eval(['Deltachi', int2str(j)]) ;
	eval(['z', int2str(j)]) ] ;
end

xinitial = [ ...
	c		
	l		
	0	
	lambda	
	lambdae	
	lambdas	
	lambdad	
	lambdai	
	Vi		
	Vd		
	Vs		
	e		
	e_by_c	
	e_by_l	
	in		
	d		
	s		
	delta 	] ;

for j=1:num_sectors
xinitial = [ xinitial ;
	eval(['c', int2str(j)]) ;
	eval(['l', int2str(j)]) ;
	eval(['m', int2str(j)]) ;
	eval(['C', int2str(j)]) ;
	eval(['L', int2str(j)]) ;
	eval(['M', int2str(j)]) ;
	eval(['chi', int2str(j), 't']) ;
	eval(['mbar', int2str(j)]) ;
	eval(['Vmbar', int2str(j)]) ;
	eval(['P', int2str(j)]) ] ;
end

% Run SIR model to get a guess of final system
[s_fin, i_fin, r_fin, d_fin] = sir(rho, gamma, delta, kappa, 0, 50, 999, 1) ;

in = i_fin(end)/1000 ;
d = d_fin(end)/1000 ;
s = s_fin(end)/1000 ;
r = r_fin(end)/1000 ;

delta = deltascale*(deltabar + exp(phi * (kappa * in))-1) ;

for j=1:num_sectors
eval(['m', int2str(j), ' = 0 ;']) ;
eval(['M', int2str(j), ' = 0 ;']) ;
eval(['l', int2str(j), ' = 1 ;']) ;
eval(['L', int2str(j), ' = (1-d-kappa*in)*l', int2str(j), ' ;']) ;
eval(['mbar', int2str(j), ' = 0 ;']) ;
eval(['chi', int2str(j), 't = chibar', ...
    int2str(j), '*(1-Deltachi', int2str(j), '*(1-exp(-rhom*mbar', int2str(j), '))) ;']) ;
eval(['M', int2str(j), ' = 0 ;']) ;
eval(['c', int2str(j), ' = (1-d-kappa*in)*(l', ...
    int2str(j), ' - chi', int2str(j), 't/2*m', int2str(j), '^2) ;']) ;
eval(['C', int2str(j), ' = (N-d)*c', int2str(j), ' ;']) ;
eval(['P', int2str(j), ' = z', int2str(j), '*(c', int2str(j), '/C', int2str(j), ')^(-1/sigma) ;']) ;
end

c = 0 ;
l = 0 ;
for j=1:num_sectors
    eval(['c = c + (1/num_sectors)*z', int2str(j), '*c', int2str(j), '^((sigma-1)/sigma) ;']) ;
    eval(['l = l + (1/num_sectors)*(1/(1+eta))*l', int2str(j), '^(1+eta) ;']) ;
end
c = c^(sigma/(sigma-1)) ;

e_by_c = 0 ;
e_by_l = 0 ;

for j=1:num_sectors
    eval(['e_by_c = e_by_c + (1/num_sectors)*(1-d)*ec', int2str(j), ...
        '*c', int2str(j), '*C', int2str(j), ';']) ;
    eval(['e_by_l = e_by_l + (1/num_sectors)*(1-d-kappa*in)*(1-m', int2str(j), ...
        ')*el', int2str(j), '*l', int2str(j), ...
        '*L', int2str(j), '*(1-M', int2str(j),') ;']) ;
end

e       = e0 + e_by_c + e_by_l ;
lambdae = ((1/c)-l1^eta)/(el1*(1-m1)*(1-M1)*L1+ec1*C1) ;
lambda  = (1/P1)*(l1^eta + lambdae*el1*(1-m1)*(1-M1)*L1) ;

lambdadtmp = 0 ;
for j=1:num_sectors
eval(['lambdadtmp = lambdadtmp - (1/num_sectors)*lambda*(l', int2str(j), ...
    '-chi', int2str(j), 't/2*m', int2str(j), '^2 - c', int2str(j), ...
    ') + lambdae*(ec', int2str(j), '*c', int2str(j), '*C', int2str(j), ...
    '+ (1-m', int2str(j), ')*el', int2str(j), '*l', int2str(j), '*L', int2str(j), '*(1-M', int2str(j), ')) ;']) ;
end

lambdad = (1/((1-(1/beta)))) * ( l^(1+eta) / (1+eta) - ud - log(c) + lambdadtmp) ;

lambdaitmp = 0 ;
for j=1:num_sectors
eval(['lambdaitmp = lambdaitmp - (1/num_sectors)*kappa*lambda*(l', int2str(j), ...
    '-chi', int2str(j), 't/2*m', int2str(j), '^2) + kappa*lambdae*((1-m', int2str(j), ')*el', int2str(j), '*l', int2str(j), '*L', int2str(j), '*(1-M', int2str(j), ')) ;']) ;
end

lambdai = (1/((1-rho-delta*kappa-1/beta))) * ( kappa * l^(1+eta) / (1+eta) - us*kappa  + lambdaitmp - delta*kappa*lambdad ) ;
lambdas = 0 ;
    
Vmtmp = 0 ;
for j=1:num_sectors
eval(['Vmtmp = Vmtmp + P', int2str(j), ...
    '*chibar', int2str(j), '*Deltachi', int2str(j), ...
    '*exp(-mbar', int2str(j), ')*m', int2str(j), '^2 ;']) ;
end
    
for j=1:num_sectors
eval(['Vmbar', int2str(j), ' = (1/(1-beta*psim))*(1/2)*(1/num_sectors)*lambda*(1-d-kappa*in)*Vmtmp ;']) ;
end

Vi = -(1/beta)*lambdai ;
Vs = -(1/beta)*lambdas ;
Vd = -(1/beta)*lambdad ;

xfinal = [ ...
	c		
	l		
	0	
	lambda	
	lambdae	
	lambdas	
	lambdad	
	lambdai	
	Vi		
	Vd		
	Vs		
	e		
	e_by_c	
	e_by_l	
	in		
	d		
	s		
	delta 	] ;

for j=1:num_sectors
xfinal = [ xfinal ;
	eval(['c', int2str(j)]) ;
	eval(['l', int2str(j)]) ;
	eval(['m', int2str(j)]) ;
	eval(['C', int2str(j)]) ;
	eval(['L', int2str(j)]) ;
	eval(['M', int2str(j)]) ;
	eval(['chi', int2str(j), 't']) ;
	eval(['mbar', int2str(j)]) ;
	eval(['Vmbar', int2str(j)]) ;
	eval(['P', int2str(j)]) ] ;
end

save input params xinitial xfinal ee_i0_init
